home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Tech Arsenal 1
/
Tech Arsenal (Arsenal Computer).ISO
/
tek-02
/
wctunits.zip
/
MATH.PAS
< prev
next >
Wrap
Pascal/Delphi Source File
|
1991-07-30
|
7KB
|
318 lines
unit math;
{ Written by William C. Thompson - wct@po.cwru.edu }
{ This unit was written to perform several basic mathematical calculations }
{ This unit automatically generates a random seed }
interface
const e=2.7182818284905;
function adjust(a:real):real;
procedure quad(a,b,c: real; var x1,x2: real; var im:real);
function relerror(observed,actual: real):real;
function pow(a,b:real):real;
function sign(r:real):integer;
function rmax(r,s:real):real;
function rmin(r,s:real):real;
function imax(m,n:integer):integer;
function imin(m,n: integer):integer;
function log(x,base: real): real;
function tan(a:real):real;
function arcsin(x:real):real;
function arccos(x:real):real;
function arctan2(y,x:real):real;
function gcd(m,n:longint):word;
function lcm(m,n:integer):integer;
function prime(n:longint):boolean;
function rnd(a,b:real):real;
function rnd2(a,b:real):real;
function gaussian(mu,sigma:real):real;
function distance(x1,y1,x2,y2:real):real;
function findline(x1,y1,x2,y2: real; var a,b:real):boolean;
function hero(a,b,c:real):real;
function factorial(n:word):extended;
function stirling(n:real):extended;
function combination(n,r:word):word;
function permutation(n,r:word):word;
implementation
function adjust(a:real):real;
{ Adjusts angle to fit into (-pi,pi] }
begin
repeat
if a<=-pi then a:=a+2*pi;
if a>pi then a:=a-2*pi
until (a>-pi) and (a<=pi);
adjust:=a
end;
procedure quad(a,b,c:real; var x1,x2: real; var im:real);
{ Solves any quadratic equation. If im=0, x1 and x2 are
two real solutions. If im<>0, x1 and x2 are real parts
of two solutions and im is imaginary part. The two
solutions would be x1 + im I and x2 - im I }
var d,q: real;
begin
im:=0.0;
if a=0 then begin
x1:=-c/b;
x2:=-c/b
end
else begin
b:=b/a;
c:=c/a;
a:=1;
d:=b*b-4*c;
q:=-b/2;
if d<0 then begin
x1:=q;
x2:=q;
im:=sqrt(-d)/2
end
else begin
if b<0 then q:=q+sqrt(d)/2
else q:=q-sqrt(d)/2;
x1:=q;
if q=0 then x2:=-b
else x2:=c/q;
if x2<x1 then begin
d:=x1;
x1:=x2;
x2:=d
end
end { d>=0 }
end { a<>0 }
end;
function relerror(observed, actual: real):real;
{ Relative error }
begin
if actual=0.0 then relerror:=abs(observed)
else relerror:=abs(observed/actual-1)
end;
function pow(a,b: real):real;
{ Computes a^b }
begin
if a=0 then
if b=0 then pow:=1 { 0^0 = 1 }
else if b<0 then pow:=exp(b*ln(a)) { force error }
else pow:=0 { 0^x = 0 }
else if a<0 then
if abs(b)<1e-10 then pow:=1
else if relerror(b,round(b))<1e-8 then
pow:=(1-2*ord(odd(round(b))))*exp(b*ln(abs(a)))
else if (relerror(1/b,round(1/b))<1e-8) and odd(round(1/b)) then
pow:=-exp(b*ln(abs(a)))
else pow:=exp(b*ln(a)) { force error }
else pow:=exp(b*ln(a))
end;
function sign(r:real):integer;
begin
if r<0 then sign:=-1
else if r>0 then sign:=1
else sign:=0
end;
function rmax(r,s:real):real;
begin if s>=r then rmax:=s else rmax:=r end;
function rmin(r,s:real):real;
begin if s<=r then rmin:=s else rmin:=r end;
function imax(m,n:integer):integer;
begin if m>=n then imax:=m else imax:=n end;
function imin(m,n:integer):integer;
begin if m<=n then imin:=m else imin:=n end;
function log(x,base:real):real;
{ Computes log of x base base }
begin
if (base=0) and (x=0) then log:=1
else if (base=0) and (x=1) then log:=0
else log:=ln(x)/ln(base)
end;
function tan(a:real):real;
begin
a:=adjust(a);
if abs(a-pi/2)<1e-10 then tan:=9.99e+37
else if abs(a-3*pi/2)<1e-10 then tan:=-9.99e+37
else tan:=sin(a)/cos(a)
end;
function arcsin(x:real):real;
begin
if x<0 then arcsin:=-arcsin(-x)
else
if x=1 then arcsin:=pi/2
else arcsin:=arctan(x/sqrt(1-x*x))
end;
function arccos(x:real):real;
begin
if x<0 then arccos:=pi-arccos(-x)
else
if x=0 then arccos:=pi/2
else arccos:=arctan(sqrt(1-x*x)/x)
end;
function arctan2(y,x:real):real;
{ Computes angle of point (x,y) }
var at2:real;
begin
if x=0 then
if y>=0 then arctan2:=pi/2
else arctan2:=-pi/2
else begin
at2:=arctan(abs(y/x));
if x>0 then
if y>0 then arctan2:=at2
else arctan2:=adjust(-at2)
else
if y>0 then arctan2:=adjust(pi-at2)
else arctan2:=adjust(-pi+at2);
end
end;
function gcd(m,n:longint):word;
{ Greatest Common Divisor }
var k,l,r:integer;
begin
m:=abs(m); n:=abs(n);
if m=0 then gcd:=n else begin
k:=m;
l:=n;
while n<>0 do begin
r:=m mod n;
m:=n;
n:=r
end;
gcd:=m
end
end;
function lcm(m,n:integer):integer;
{ Least common multiple }
begin
if (m=0) and (n=0) then lcm:=0
else lcm:=abs(round((m*1.0)*n) div gcd(m,n))
end;
function prime(n:longint):boolean;
{ Tests a number for primeness }
var
p: boolean;
i, limit: word;
begin
n:=abs(n);
if n in [0..1] then prime:=false
else if n in [2..3] then prime:=true
else begin
p:=false;
i:=2;
limit:=round(sqrt(n));
while (i<=limit) and not p do begin
if (n mod i=0) then p:=true;
i:=i+1
end;
prime:=not p
end
end;
function rnd(a,b:real):real;
{ divides interval from a to b, inclusive, into 65535 values and
chooses one }
begin
rnd:=a+random(65535)/65534*(b-a)
end;
function rnd2(a,b: real):real;
{ Similar to rnd, but MANY more divisions (over 4 million) }
begin
rnd2:=a+(65535.0*random(65535)+random(65535))/(sqr(65535.0)-1)*(b-a)
end;
function gaussian(mu,sigma:real):real;
{ Produces a Gaussian distributed random variable with mean mu
and standard deviation sigma }
var
r,v1,v2: real;
begin
repeat
v1:=rnd(-1,1);
v2:=rnd(-1,1);
r:=sqr(v1)+sqr(v2)
until (r<1) and (r>0);
gaussian:=v1*sqrt(-2*ln(r)/r)*sigma+mu
end;
function distance(x1,y1,x2,y2:real):real;
begin
distance:=sqrt(sqr(x1-x2)+sqr(y1-y2))
end;
function findline(x1,y1,x2,y2: real; var a,b:real):boolean;
{ Finds equation for line y=ax+b between (x1,y1) and (x2,y2) }
begin
if x2=x1 then begin
findline:=false;
exit
end;
findline:=true;
a:=(y2-y1)/(x2-x1);
b:=(x2*y1-x1*y2)/(x2-x1);
end;
function hero(a,b,c:real):real;
{ Finds area of triangle with sides of length s }
var
s: real;
begin
s:=a/2+b/2+c/2;
hero:=sqrt(s*(s-a)*(s-b)*(s-c));
end;
function factorial(n:word):extended;
{ computes factorial }
var
i: word;
f: extended;
begin
f:=1;
for i:=1 to n do f:=f*i;
factorial:=f
end;
function stirling(n:real):extended;
{ computes factorial according to Stirling's approximation -
this may not be correct... it should be, but ask your math professor }
var
s: extended;
begin
s:=1/n;
s:=s*exp(n*ln(n));
stirling:=s*exp(-n)/sqrt(2*pi);
end;
function combination(n,r:word):word;
{ computes nCr }
begin
combination:=round((factorial(n)/factorial(n-r))/factorial(r))
end;
function permutation(n,r:word):word;
{ computes nPr }
begin
permutation:=round(factorial(n)/factorial(n-r))
end;
begin
randomize
end.